val = zeros(nSamp,numel(lambda));
val(:,ii) = iePoisson(lambda(ii),'nSamp',[nSamp,1]);
h = histogram(val(:,1),'Normalization','pdf'); hold on;
h.BinLimits = [-5 40]; h.BinWidth = 1; % h.FaceColor = [0.8 0.8 0.8];
h = histogram(val(:,2),'Normalization','pdf');
h.BinLimits = [-5 40]; h.BinWidth = 1; % h.FaceColor = [0.5 0.5 0.5];
h = histogram(val(:,3),'Normalization','pdf');
h.BinLimits = [-5 40]; h.BinWidth = 1; % h.FaceColor = [0.2 0.2 0.2];
ylabel('Probability density')
val = zeros(nSamp,numel(lambda));
val(:,ii) = sqrt(lambda(ii))*randn(nSamp,1) + lambda(ii);
h = histogram(val(:,1),'Normalization','pdf'); hold on;
h.BinLimits = [-5 40]; h.BinWidth = 1; % h.FaceColor = [0.8 0.8 0.8];
h = histogram(val(:,2),'Normalization','pdf');
h.BinLimits = [-5 40]; h.BinWidth = 1; % h.FaceColor = [0.5 0.5 0.5];
h = histogram(val(:,3),'Normalization','pdf');
h.BinLimits = [-5 40]; h.BinWidth = 1; % h.FaceColor = [0.2 0.2 0.2];
ylabel('Probability density')
theMosaic = mosaicLoad('cmosaic_0.5-0.5_0.0-0.0');theMosaic.wave = 400:1:700;
theLens.wave = 400:1:700;
theMosaic.macular.density = 0.35; % 0.35 is Default
theMosaic.macular.eccDensity(0:1:15);
theLens.density = 1; % 1 is default
ieNewGraphWin; plot(wave,theMosaic.pigment.quantaFundamentals);
ieNewGraphWin; plot(wave,theMosaic.qe);
ieNewGraphWin([],'wide');
theMosaic.macular.density = theMosaic.macular.eccDensity(0);
coneFundamentals = bsxfun(@times, theLens.transmittance', theMosaic.qe);
plot(wave,coneFundamentals,'Linewidth',3);
xlabel('Wavelength (nm)');
legend('L-cones','M-cones','S-cones');
theMosaic.macular.density = theMosaic.macular.eccDensity(10); % 0.35 is Default
coneFundamentals = bsxfun(@times, theLens.transmittance', theMosaic.qe);
plot(wave,coneFundamentals,'Linewidth',3);
xlabel('Wavelength (nm)');
title('10 deg eccentricity');
% I stopped the ISETBioCSFGenerator inside the
% cMosaic.oiEnsembleGenerate
% method. While there, I saved the key optics variables in the 'projects'
% save('oiData','obj','oiSamplingGridDegs','varargin');
load('oiData','obj','oiSamplingGridDegs','varargin');
% See oiPosition, not this!
% The key variable is obj, the ISETBio cMosaic
% obj.oiEnsembleGenerate runs to create the psfs given the parameters
% edit cMosaic.oiEnsembleGenerate
% The meaning of the parameters is:
% centerPos = oiSamplingGridDegs;
{'zernikeDataBase'} {'Polans2015'}
{'pupilDiameterMM'} {[3]}
{'subtractCentral…'} {[0]}
{'flipPSFUpsideDown'} {[1]}
{'wavefrontSpatia…'} {[301]}
% Convert the varargin to a struct that we can edit
for ii=1:2:numel(varargin)
args.(varargin{ii}) = varargin{ii+1};
mp(1:4,:) = repmat([0.5 0.5 0.5],[4,1]);
centers = [2,0; 8,0; 20,0];
ieNewGraphWin([],'wide');
centerPos = centers(ii,:);
args.zernikeDataBase = 'Artal2012'; % args.zernikeDataBase = 'Polans2015';
args.pupilDiameterMM = 4;
args.subtractCentralRefraction = false;
args.zeroCenterPSF = true;
args.flipPSFUpsideDown = true;
args.wavefrontSpatialSamples = 301;
[oiEnsemble, psfEnsemble] = obj.oiEnsembleGenerate(centerPos,args);
thisPSF = psfEnsemble{1};
thisW = theW(jj); % Select a wavelength, % 3 is 450 nm, 7 is 550 nm
% h = mesh(thisPSF.supportX,thisPSF.supportY,thisPSF.data(:,:,thisW)/max2(thisPSF.data(:,:,thisW)));
imagesc(thisPSF.supportX,thisPSF.supportY,thisPSF.data(:,:,thisW)/max2(thisPSF.data(:,:,thisW)));
axis image; colormap(mp);
tMarks = (-15:5:15); mn = -8; mx = 8;
set(gca,'xtick',tMarks,'xlim',[mn mx],'ylim',[mn mx]);
zlabel('Relative intensity');
% title(sprintf('Wave %d Sub %d',thisPSF.supportWavelength(thisW),args.subjectID));
scene = sceneCreate('lstar');
scene = sceneSet(scene,'fov',5);
oi = oiCompute(oi,scene);
cmP.spatialDensity = cDensity;
cmP.integrationTime = eTime; % 10 ms
photonNoise = cmP.absorptions;
mx = max(photonNoise(:));
ieNewGraphWin; imagesc(photonNoise,[0 mx]);
axis image; colormap(hot(128));
title('Poisson');brighten(0.3);
thisRow = photonNoise(400,:);
p = plot(thisRow,'LineWidth',1,'Color',[1 1 1]);
ylabel('Excitations'); xlabel('Position'); grid on;
cmG.spatialDensity = cDensity;
cmG.integrationTime = eTime;
noNoise = cmG.absorptions;
equivNoise = randn(size(noNoise))*sqrt(mn);
gaussNoise = noNoise + equivNoise;
gaussNoise(gaussNoise<0) = 0;
ieNewGraphWin; imagesc(gaussNoise,[0 mx]);
axis image; colormap(hot(128));
title('Gaussian'); brighten(0.3);
thisRow = gaussNoise(400,:);
p = plot(thisRow,'LineWidth',1,'Color',[1 1 1]);
ylabel('Excitations'); xlabel('Position'); grid on;
cmG.absorptions = gaussNoise;